\documentclass[english,uzh-en]{article}

\usepackage[T1]{fontenc}
\usepackage[latin9]{inputenc}
\usepackage{babel}
\usepackage{graphicx, lipsum}
\usepackage[a4paper, total={17cm, 25cm}]{geometry}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{hyperref}

\begin{document}

\vspace*{1cm}
\noindent {\LARGE \textbf{Longitudinal Study CD8-Tcell Report}}\\

<<echo=FALSE, warning=FALSE>>=
rm(list = ls())

library(readxl, quietly = TRUE, warn.conflicts = FALSE)
library(xtable, quietly = TRUE, warn.conflicts = FALSE)
# install.packages("biostatUZH", repos="http://R-Forge.R-project.org")
library(biostatUZH, quietly = TRUE, warn.conflicts = FALSE)
library(lme4, quietly = TRUE, warn.conflicts = FALSE)
library(lmerTest, quietly = TRUE, warn.conflicts = FALSE)



PiB <- as.data.frame(read_excel("../Data/PiB Imaging Data Kirabali-Gericke.xlsx"))
BCell_data <- as.data.frame(read_excel("../Data/B cells_longitudinal data_re-arranged_Nov2020.xlsx"))
CD4TCell_data <- as.data.frame(read_excel("../Data/CD4_longitudinal T cell data_CG_2020_03_20.xlsx"))
CD8TCell_data <- as.data.frame(read_excel("../Data/CD8_longitudinal T cell data_CG_2020_03_20.xlsx"))


source("eda.R")

### PiBGroup as in BCell_data
getoriginalPIB <- function(subject = "01_4710_01") {
  out <- BCell_data[BCell_data$Subject == subject, "PiBGroup"]

  if (length(unique(out)) == 1)
    return(unique(out))
  else
    stop("not unique groupallocation")
}

noteqsub <- NULL
for (i in 1:dim(CD8TCell_data)[1]) {
  group <- getoriginalPIB(subject = as.character(CD8TCell_data$Subject[i]))

  if(CD8TCell_data$PiBGroup[i] != group) {
    noteqsub <- c(noteqsub, as.character(CD8TCell_data$Subject[i]))
  }

  CD8TCell_data$PiBGroup[i] <- group
}

### PiBGroup as in BCell_data
getoriginalPIBvalues <- function(subject = "01_4710_01", timep = 1) {
  out <- BCell_data2[BCell_data2$Subject == subject & BCell_data2$Timepoint == timep, "PiB"]

  if (length(unique(out)) == 1)
    return(unique(out))
  else
    stop("not unique groupallocation")
}

noteqsubpibvalue <- NULL
for(j in 1:3) {
  for (i in 1:dim(CD8TCell_data)[1]) {
    pib <- getoriginalPIBvalues(subject = as.character(CD8TCell_data$Subject[i]), timep = j)

    if (j !=2) {
      if(unique(CD8TCell_data[CD8TCell_data$Subject ==as.character(CD8TCell_data$Subject[i]) & CD8TCell_data$Timepoint == j, "PiB"]) != pib) {
        noteqsubpibvalue <- c(noteqsubpibvalue, as.character(unique(CD8TCell_data$Subject[i])))
      }
    }

    CD8TCell_data[CD8TCell_data$Subject ==as.character(CD8TCell_data$Subject[i]) & CD8TCell_data$Timepoint == j, "PiB"] <- pib

  }
}
@


<<cache=FALSE, echo=FALSE, warning=FALSE>>=
expit <- function(input) {
  exp(input)/(1 + exp(input))
}

library(lme4)

getCell2Cluster <- function(id) {
  clusterdata <- NULL
  for (i in id) {
    clusterdata <- rbind(clusterdata, CD8TCell_data[CD8TCell_data$Cluster_id == i, ])
  }

  return(clusterdata)
}

Clustersubsetfor2t <- function(clusterid = i, grouptwo = TRUE) {
  CD8TCell_cli <- getCell2Cluster(id = clusterid)
  if (grouptwo) {
      CD8TCell_cli <- CD8TCell_cli[CD8TCell_cli$PiBGroup == 2, ]
  }

  pibsubjects <- unique(CD8TCell_cli$Subject)
  subsetfor2 <- function(subj = pibsubjects[1]) {
    subjdata <- CD8TCell_cli[CD8TCell_cli$Subject == subj, ]

    if (!is.na(subjdata[subjdata$Timepoint == 1, "Abundance"]) & !is.na(subjdata[subjdata$Timepoint == 3, "Abundance"])) {
      return((subjdata[-which(subjdata$Timepoint == 2), ])) }

    if (is.na(subjdata[subjdata$Timepoint == 1, "Abundance"]) & !is.na(subjdata[subjdata$Timepoint == 3, "Abundance"])) {
      return((subjdata[-which(subjdata$Timepoint == 1), ])) }

    if (!is.na(subjdata[subjdata$Timepoint == 1, "Abundance"]) & is.na(subjdata[subjdata$Timepoint == 3, "Abundance"])) {
      return((subjdata[-which(subjdata$Timepoint == 3), ])) }

    return((subjdata))
  }

  cli2time_data <- CD8TCell_cli[-c(1:dim(CD8TCell_cli)[1]), ]
  for (i in 1:length(pibsubjects)) {
    cli2time_data <- rbind(cli2time_data, subsetfor2(subj = pibsubjects[i]))
  }

  return(cli2time_data)
}

simple_normalfit <- function(clusterid = 1) {
  CD8TCell_sm_data <- Clustersubsetfor2t(clusterid = clusterid, grouptwo = F)

  tmpsubj <- unique(CD8TCell_sm_data$Subject)
  thisDelta_Ab <- sapply(tmpsubj, function(x) {
    tmpsubdata <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, c("Timepoint", "Abundance")]
    if (dim(tmpsubdata)[1] < 2)
      return(NA)

    return(tmpsubdata$Abundance[2] - tmpsubdata$Abundance[1])
  })

  thisAge <- sapply(tmpsubj, function(x) {
    tmpsubjage <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, "Age"]
    if(!is.na(tmpsubjage[1])) {
       return(tmpsubjage[1])
    } else {
      return(NA) }
    })

  thisPiBGroup <- sapply(tmpsubj, function(x) {
    tmpsubjgroup <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, "PiBGroup"]
    return(tmpsubjgroup[1]) })

  thisDelta_Ab <- scale(log(thisDelta_Ab + sqrt(thisDelta_Ab^2 + 1)))

  sm_fit <- lm(thisDelta_Ab ~ 1 + thisPiBGroup + thisAge)

  return(sm_fit)
}

all_simple_normalfits <- list(NULL)
for (i in 1:21) {
  all_simple_normalfits[[i]] <- simple_normalfit(clusterid = i)
}

simple_binomialfit <- function(clusterid = 1) {
  CD8TCell_sm_data <- Clustersubsetfor2t(clusterid = clusterid, grouptwo = F)

  tmpsubj <- unique(CD8TCell_sm_data$Subject)
  thisDelta_Ab <- sapply(tmpsubj, function(x) {
    tmpsubdata <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, c("Timepoint", "Abundance")]
    if (dim(tmpsubdata)[1] < 2)
      return(NA)

    return(tmpsubdata$Abundance[2] - tmpsubdata$Abundance[1])
  })

  thisAge <- sapply(tmpsubj, function(x) {
    tmpsubjage <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, "Age"]
    if(!is.na(tmpsubjage[1])) {
       return(tmpsubjage[1])
    } else {
      return(NA) }
    })

  thisPiBGroup <- sapply(tmpsubj, function(x) {
    tmpsubjgroup <- CD8TCell_sm_data[CD8TCell_sm_data$Subject == x, "PiBGroup"]
    return(tmpsubjgroup[1]) })

  thisDelta_Ab <- scale(log(thisDelta_Ab + sqrt(thisDelta_Ab^2 + 1)))

  sm_fit <- glm(thisPiBGroup ~ 1 + thisDelta_Ab + thisAge, family = binomial(link = "logit"))

  return(sm_fit)
}

all_simple_binomfits <- list(NULL)
for (i in 1:21) {
  all_simple_binomfits[[i]] <- simple_binomialfit(clusterid = i)
}

mixed_normalfit2timepoints <- function(clusterid = 11) {

  cli2time_data <- Clustersubsetfor2t(clusterid = clusterid, grouptwo = F)

  mm2_fit <- lmerTest::lmer((PiB) ~ 1 + (Abundance*PiBGroup) + (Age) + (1|Subject), data = cli2time_data)

  return(mm2_fit)
}

all_mixed_normalfit2timepoints <- list(NULL)
for (i in 1:21) {
  all_mixed_normalfit2timepoints[[i]] <- mixed_normalfit2timepoints(clusterid = i)
}

trans_mixed_normalfit2timepoints <- function(clusterid = 11) {

  cli2time_data <- Clustersubsetfor2t(clusterid = clusterid, grouptwo = F)
  cli2time_data$Abundance <- scale(log(cli2time_data$Abundance + sqrt(cli2time_data$Abundance^2 + 1)))

  mm2_fit <- lmerTest::lmer((PiB) ~ 1 + (Abundance*PiBGroup) + (Age) + (1|Subject), data = cli2time_data)

  return(mm2_fit)
}

all_trans_mixed_normalfit2timepoints <- list(NULL)
for (i in 1:21) {
  all_trans_mixed_normalfit2timepoints[[i]] <- trans_mixed_normalfit2timepoints(clusterid = i)
}

@


<<echo=FALSE, fig.height=7, eval=FALSE>>=
{
  par(mfrow = c(4, 1))
  plot(sapply(1:21, function(i) { AIC(all_simple_binomfits[[i]]) }), type = "b", col = "blue", xlab = "cluster id", ylab = "AIC", main = "simple binom")
  plot(sapply(1:21, function(i) { AIC(all_simple_normalfits[[i]]) }), type = "b", col = "blue", xlab = "cluster id", ylab = "AIC", main = "simple normal")
  plot(sapply(1:21, function(i) { AIC(all_mixed_normalfit2timepoints[[i]]) }), type = "b", col = "blue", xlab = "cluster id", ylab = "AIC", main = "mixed 2 timepoints normal")
  plot(sapply(1:21, function(i) { AIC(all_trans_mixed_normalfit2timepoints[[i]]) }), type = "b", col = "blue", xlab = "cluster id", ylab = "AIC", main = "mixed transformed 2 timepoints normal")
}
@

<<echo=FALSE,warning=FALSE,message=FALSE,results='asis', fig.height=4.5, cache=FALSE>>=

for(i in 1:21) {
  if (FALSE) {
    tmptable <- biostatUZH::tableRegression(all_simple_binomfits[[i]], xtable = FALSE)
    tmptable <-  cbind(tmptable, matrix(summary(all_simple_binomfits[[i]])$coefficients[2:3,2], ncol = 1))
    colnames(tmptable)[4] <- "se"
    print(xtable::xtable(tmptable, caption = paste("Simple binomial model for cluster", i)))

    { par(mfrow = c(1,3))
      for (j in 1:2)
        plot(all_simple_binomfits[[i]], which = j)

      x <- seq(range(all_simple_binomfits[[i]]$data$thisDelta_Ab, na.rm = T)[1],
               range(all_simple_binomfits[[i]]$data$thisDelta_Ab, na.rm = T)[2], length.out = 100)
      plot(x = x,
           y = expit(coef(all_simple_binomfits[[i]])[1] + coef(all_simple_binomfits[[i]])[2]*x +
                       coef(all_simple_binomfits[[i]])[3]*mean(all_simple_binomfits[[i]]$data$thisAge)),
           type = "l", xlab = "Delta Abundance", ylab = "Prob for Group 2",
           ylim = c(0,1))
      lines(x = x, y = expit(coef(all_simple_binomfits[[i]])[1] + coef(all_simple_binomfits[[i]])[2]*x +
                               coef(all_simple_binomfits[[i]])[3]*min(all_simple_binomfits[[i]]$data$thisAge)), col = "red")
      lines(x = x, y = expit(coef(all_simple_binomfits[[i]])[1] + coef(all_simple_binomfits[[i]])[2]*x +
                               coef(all_simple_binomfits[[i]])[3]*max(all_simple_binomfits[[i]]$data$thisAge)), col = "green")
      legend("topleft", c("mean age", "min age", "max age"), col = c("black", "red", "green"), pch = 16)
    }

    tmpregtable <- biostatUZH::tableRegression(all_simple_normalfits[[i]], xtable = FALSE)
    tmpregtable <-  cbind(tmpregtable, matrix(summary(all_simple_normalfits[[i]])$coefficients[,2], ncol = 1))
    colnames(tmpregtable)[4] <- "se"
    print(xtable::xtable(tmpregtable, caption = paste("Simple binomial model for cluster", i)))

    { par(mfrow = c(1,3))
      for (j in 1:2)
        plot(all_simple_normalfits[[i]], which = j)
    }
  }

  if (FALSE) {
    modelfitsoutput(listoffits = all_mixed_normalfit2timepoints[[i]], Clustersubsetfor2t(i), clusterpathimputed = T,
                    tablecaption = "mixed effects 2timepoints") }

    modelfitsoutput(listoffits = all_trans_mixed_normalfit2timepoints[[i]], Clustersubsetfor2t(i), clusterpathimputed = T,
                    tablecaption = "mixed effects transformed 2timepoints")
  cat("\\newpage")

}

@

\end{document}
